Multi-particle dispersion in fully developed turbulence 
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The statistical geometry of dispersing Lagrangian clusters of four particles (tetrahedra) is stud- 
ied by means of high-resolution direct numerical simulations of three-dimensional homogeneous 
isotropic turbulence. We give the first evidence of a self-similar regime of shape dynamics char- 
acterized by almost two-dimensional, strongly elongated geometries. The analysis of four-point 
velocity-difference statistics and orientation shows that inertial-range eddies typically generate a 
straining field with a strong extensional component aligned with the elongation direction and weak 
extensional/compressional components in the orthogonal plane. 

PACS numbers: 47.27-i 



One of the most characteristic attributes of turbulence 
is the efficient dispersion and mixing of advected La- 
grangian particles [1] . Even though turbulent dispersion 
bears some similarities to Brownian motion, especially 
at very large scales and for long times, it has a much 
richer structure at small scales. This is already visible 
at the level of single particle dispersion, which is char- 
acterized by non-trivial time-correlations of the velocity 
experienced by the particle along its trajectory (see e.g. 
[2]). The statistics of pair dispersion display interest- 
ing properties as well (see e.g. [3, 4, 5]), yet the com- 
plexity of Lagrangian turbulence is particularly evident 
when looking at the dispersion of three or more particles. 
This calls for the description of the geometrical properties 
of Lagrangian dispersion - the "shape" of the particles' 
cloud as well as its "size". The geometrical characteri- 
zation of dispersion proved extremely important for the 
understanding of the problem of passive scalar advec- 
tion [6] and provides the basis for the efficient modelling 
of the small-scale velocity dynamics itself [7]. Previous 
studies dealt with two-dimensional flows [8, 9], synthetic 
flows [10, 11] or three-dimensional turbulence at moder- 
ate Reynolds numbers [7, 12, 13]. In this Letter we study 
multi-particle Lagrangian statistics by means of high res- 
olution direct numerical simulations of three-dimensional 
Navier-Stokes turbulence. Simulations were done at res- 
olutions of 1024 3 corresponding to a Reynolds number 
R\ ~ 280 (see Ref. [14]). The other parameters of the 
numerical simulation are as follows: energy dissipation 
e = 0.81(8), viscosity v — 8.8 • 10~ 4 , Kolmogorov length 
scale 77 = 5 • 10~ 3 , integral scale L — 3.14, Lagrangian 
velocity autocorrelation time Tl = 1.2, Kolmogorov time 
scale t v = 3.3 • 10~ 2 . 

With the present choice of parameters the dissipa- 
tive range of length scales is exceptionally well resolved. 
Upon having reached a statistically stationary velocity 



field, the Lagrangian tracers were seeded in the flow. 
Their trajectories were integrated according to 



dx 
~dt 



u(x(t),i) 



over a time lapse of the order of a few Lagrangian corre- 
lation times, Tl. The velocity field, u, results from the 
time-integration of the three-dimensional Navier-Stokes 
equations (for further details see Ref. [14]). 

A set of 3.84 • 10 5 particles were initially seeded in 
quadruplets forming 9.6 • 10 4 regular tetrahedra of the 
size of the Kolmogorov scale, with centers of mass uni- 
formly distributed over the domain. The evolution of 
the separations between different particles in each tetra- 
hedron provides a way to quantify the shape evolution. 
As particles move with the flow the size of the tetra- 
hedra grows in time and their shape deforms, generat- 
ing a variety of irregular objects. A description of this 
process is then given in terms of the probability density 
functions (pdf) of sizes and shapes. Within the inertial 
range of scales a self-similar evolution of size according 
to Richardson's law and a stationary shape distribution 
are expected. Figure 1 shows a sample of the tetrahe- 
dra evolving in the turbulent flow. The presence of very 
different shapes, from almost regular to very flat and 
elongated involving the interaction of diverse scales, is 
evident. 

In order to characterize the shape dynamics quanti- 
tatively, it is useful to introduce the following change 
of coordinates [7]: po = (xi + X2 + xj, + Xi)/2 pi = 
(x 2 - Xi)/y/2, pi = {2x 3 - x 2 - xx)/VH, P3 = (3^4 - 
x 3 — x 2 — x\)/ y/12. By virtue of the statistical homogene- 
ity of the velocity field as well as of the initial distribution 
of the centers of mass, the Lagrangian statistics do not 
depend on po. The information about the particle sep- 
arations can be embodied in the square matrix p whose 
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FIG. 1: Snapshot at t = 0.65T L of 480 tetrahedra evolving 
in the turbulent flow starting from regular tetrahedra at the 
Kolmogorov scale. 



columns are the three vectors pi with i = 1, 2, 3. Denot- 
ing by <?i (<7i > 32 > 33) the eigenvalues of the moment 
of inertia matrix, I = pp T (that is positive defined), 
we have that the size of the tetrahedron is r = y/tr{I) = 

V.9i + 32 + 93 = \JlJ2i,j \ x i ~ x j\ 2 > whereas the volume 
can be expressed as V — | det(p) = |y '313233- A con- 
venient characterization of shapes is given in terms of 
the dimcnsionlcss quantities Ii = gi/r 2 (where obviously 
h + h + h = 1). For a regular tetrahedron one has 
h = I2 = I3 — 1/3. If the four points are coplanar one 
has I3 = and for a collinear configuration I2 = I3 = 0. 

Figure 2 shows the temporal evolution of the mean 
eigenvalues of pp T for the smallest regular tetrahedra 
with 3i(0) = Sx 2 /2. Two very different regimes are evi- 
dent: at small times t < the evolution of tetrahedra is 
governed by the dissipative range of turbulence. Because 
of the smoothness and incompressibility of the velocity 
field in this range, the volume of each tetrahedron is ap- 
proximately preserved and so is its average value which 
is shown in Fig. 2. In the viscous range the shape dy- 
namics are essentially characterized by the Lagrangian 
Lyapunov exponents [15]: as a consequence the mean 
square separation r 2 grows exponentially in time. From 
the average growth rate of the logarithms of the separa- 
tions, R(t) — \pi\, areas A(t) — ^\pi X pz\ an d volumes 
V(t) = i \pi x P2 ■ p3 1 at small times, we can obtain an es- 
timation of the Lagrangian Lyapunov spectrum as shown 
in Fig. 2. We found two positive Lyapunov exponents, 
with Ait,, ~ 0.12 and A2 — Ai/4, in agreement with pre- 
vious findings at lower R\ [16, 17]. The sum of the three 
Lyapunov exponents so obtained is close to zero for times 

Up tO 3Trf. 

The exponential growth brings particle separations 
outside the dissipative range, where the velocity field be- 
comes rough and the inertial range sets in. According to 
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FIG. 2: Evolution of the mean eigenvalues g\ (+), gi (x) 
and (?3 (*) of the moment of inertia matrix i" = pp T . The 
line represents the dimensional scaling t 3 . In the inset, from 
top to bottom: evolution at small times of (In A(t)} (surface), 
(In R(t)) (distance), (lnV(t)) (volume). The linear slopes of 
the three curves in the range of times t v < t < 3t v yield 
Ai + A2, Ai and Ai + A2 + A3, respectively. 



the Kolmogorov-Richardson scaling, eigenvalues should 
grow as gi ~ i 3 . As previously reported [7], it is hard 
to extract a clear scaling regime for the shape dynamics 
shown in Fig. 2. The main reason for the lack of self- 
similarity is due to the contamination of the inertial range 
by the dissipative range. Indeed, because of the strong 
shape distortion taking place at the crossover between 
the dissipative and inertial ranges (as shown in Fig. 2 by 
the separation of the three eigenvalues), a significant frac- 
tion of tetrahedra has one side in the dissipative range 
even at times much larger than t v . In order to overcome 
this problem we have utilized the technique of doubling 
time statistics that has already been succesfully used to 
remove contaminations in the statistics of pair dispersion 
[4, 5, 18]. Here, we focus on the doubling times of the 
eigenvalues gc we compute the times, T(g{), taken by 
a tetrahedron to increase its value of gi by a factor a. 
The result is shown in Fig. 3. The presence of a scaling 
range T ~ 3 1 / 3 is more clear and the self-similarity is 
made evident by superimposing the three curves on top 
of each other by a simple multiplicative factor on the g- 
axis. The ratio of the three eigenvalues in the scaling 
range is 31 : 32 : 33 = 40 : 8 : 1, corresponding to shape 
indices I2 ~ 0.16 and I3 w 0.02. The presence of a range 
where the doubling times for different eigenvalues are the 
same is equivalent to stating that the typical shape of the 
tetrahedron is preserved while its size increases according 
to Richardson's law. 

In view of the existence of a self-similar regime for 
shape evolution, one would expect that the statistics of 
the shape indices, Ii, should reach a time- independent 
distribution. However, a direct inspection of the data 
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FIG. 3: Doubling times for the eigenvalues, (?;, of the moment 
of inertia matrix, pp T . In the inset: the same data rescaled 
on the horizontal axis with the proportions g\ : g2 : <73 = 40 : 
8:1 
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FIG. 4: Probability density function of shape indices I2 and 
73 (inset) at times t = 35t v (+) and t — 63r, ; (x). The 
full lines are the pdfs for independent, Gaussian distributed 
particle positions. 



does not support this conclusion (not shown here, the 
results do not present an appreciable scaling range in 
time in spite of the relatively high R\ as compared with 
Ref. [7]). Once more this lack of a scaling range in the 
time domain can be traced back to the contamination by 
the dissipative range dynamics. 

This difficulty can be overcome by selecting those 
tetrahedra with eigenvalues in the ranges 5- 10 2 f? 2 < g% < 
5 • 10V, 5 • 10 V < g 2 < 5 • 10V, 5 • rj 2 < g 3 < 5 • 10V ■ 
The thresholds are obtained by identifying the scaling 
ranges in Fig. 3. This procedure removes about 60% of 
the initial tetrahedra, mostly because 33 falls below its 
lower threshold. The probability density functions of the 
shape indices I2, 13 after the selection are shown in Fig. 4. 
The existence of an invariant regime appears now very 
clearly. In this regime, the normalised probability density 
functions at different times collapse, and the mean values 
hence display a plateau in time: for the third index, the 
mean value (I3) ~ 0.011 ± 0.001 is not too far from the 
Gaussian value 0.03, while the second index is concen- 
trated on values much smaller ({I2} — 0.135 ± 0.003 as 
opposed to 0.22). Those values indicate a relative abun- 
dance of flat and elongated configurations. The tendency 
to form almost two-dimensional structures has mostly 
an "entropic" origin: indeed there is a large number of 
pancake-like tetrahedra (very small I3) already for Gaus- 
sian, independent particle positions, as shown by the cor- 
responding distribution in Fig. 4. However, it has to be 
remarked that the pdf of I3 is significantly more peaked 
at small values than the Gaussian one. The preference 
for elongated structures (I2 <C Ii) has a clear dynamical 
origin, since it has no equivalent in the Gaussian ensem- 
ble. 

An interesting issue that we do not address here is 
connected with the possibility of subleading, anomalous 



scaling in the tetrahedra distribution. In the simpler case 
of particles advected by a Gaussian and white- in-time ve- 
locity field, it is known that the asymptotic behaviour of 
the multi-particle pdf, when the intial points are close, is 
governed by an expansion in zero modes and slow modes 
of a given evolution operator [6]. There, anomalous cor- 
rections emerge as sub-leading terms to the Richardson 
scaling. These corrections are connected to the anoma- 
lous scaling of the structure functions of a passive scalar 
field advected by the flow. Here, in the presence of a real 
turbulent flow, one can only argue that similar properties 
may still hold [8]. In order to check this, one should per- 
form a delicate compensation between the evolution of 
the pdf with different initial tetraehdra shapes, in order 
to cancel the leading scaling terms and to highlight the 
sub-leading contributions. 

The dynamics of the shape evolution can be eluci- 
dated by analyzing the local geometrical properties of 
Lagrangian velocities. In analogy with the relative co- 
ordinates p, we introduce the relative velocity matrix 
W: W x = {u 2 - W 2 = (2m 3 - U2 - ui)/V6, 

W3 = (3it4 — it 3 — U2 — iti)/VT2. Obviously, p = W . 
The geometrical aspects of Lagrangian velocity evolution 
can be described by the tetrahedron "turbulent diffusion" 
tensor 

K^ l -±pp T = l -{Wp T + P W T ). (1) 

The trace tr(K) = | Yl % j ( u i ~ u j) ' ( x i — x j) i s propor- 
tional to the longitudinal velocity difference multiplied by 
the separation averaged over all pairs within the tetrahe- 
dron. The geometrical information about the Lagrangian 
velocity fluctuations may be obtained from the eigenval- 
ues n\ > K2 > K3 of K which are shown in Fig. 5. On 
dimensional grounds these should grow in time as t 2 or, 
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flows. The joint effort on the numerical and experimen- 
tal side can shed further light on the geometrical statis- 
tics of Lagrangian turbulence. This, in turn, will lead 
to the development of new, more effective parameteriza- 
tions of small-scale turbulence, a problem of paramount 
importance for geophysical and industrial applications. 

We are grateful to A. Pumir for discussions and use- 
ful suggestions. The simulations were performed on the 
IBM-SP4 of the Cineca (Bologna, Italy) . We are grateful 
to C. Cavazzoni and G. Erbacci for resource allocation 
and precious technical assistance. 



FIG. 5: Evolution of the mean eigenvalues of the "turbulent 
diffusion tensor", K, as a function of the tetrahedron size, r. 
In the inset, the eigenvalues as a function of time. 

equivalently, with the tetrahedron size, r, as r 4 / 3 : this is 
satisfied to a good accuracy for all three eigenvalues, es- 
pecially as a function of size. The third eigenvalue, K3, is 
negative (notice that, strictly speaking, this makes abu- 
sive the definition of K as a diffusion tensor): geometri- 
cally this means that the local velocity field experienced 
by the tetrahedron has two extensional components, a 
strong one and a weak one, Ki S> k%, with the latter 
smaller by a factor of ten than the former, and a weak 
compressional component ~ K2- It is also interest- 
ing to study the relative orientation of the eigenvectors 
of the matrix I — pp T , i.e. the principal axes of inertia, 
and the eigenvectors of the matrix K. We found that 
the directions of the eigenvectors associated with g\ and 
Ki are preferentially aligned. About 45% of the tetrahe- 
dra show a relative angle smaller than 7r/6 (for a uniform 
distribution on a unit sphere one would have 13%). This 
agrees with the intuitive idea that strongly extensional 
velocity differences result in intense elongations approx- 
imately in the same direction. In the plane orthogonal 
to the first principal axis of inertia, the eigenvectors of I 
and K associated with the smaller eigenvalues are also 
aligned albeit to a lesser degree (about 25% of relative 
angles below 7r/6). 

The overall geometrical picture that emerges is the fol- 
lowing: tetrahedra tend to be elongated, almost copla- 
nar objects, subject to a straining velocity field that 
has a strong extensional part in the direction of elon- 
gation and relatively weak compressive and extensional 
contributions in the orthogonal plane of approximately 
equal magnitude. The recent advances in experimental 
techniques for particle tracking should soon allow pre- 
cise measurements of shape dynamics in real turbulent 
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